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Abstract - We present a new type of the EnKF for data 
assimilation in spatial models that uses diagonal approxima- 
tion of the state covariance in the wavelet space to achieve 
adaptive localization. The efficiency of the new method is 
demonstrated on an example. 
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I. Introduction 

Incorporating new data into computations in progress is 
a well-known problem in many areas, including weather 
forecasting, signal processing, and computer vision. Se- 
quential Bayesian techniques based on the state- space 
model are known as filtering or data assimilation (T), 
(2). The probability distribution of the system state is 
advanced in time by the computational model, while 
data is incorporated from time to time by modifying 
the probability distribution of the state by an application 
the Bayes theorem. Gaussian probability distributions are 
represented by their mean and covariance. An assumed 
constant state covariance then yields the classical optimal 
statistical interpolation (OSI). The Kalman filter (KF) 
evolves the state covariance, but it needs to maintain 
the covariance matrix, so it is not suitable for high- 
dimensional systems. The ensemble Kalman filter (EnKF) 
(2) replaces the state covariance by sample covariance 
of an ensemble of simulations. The EnKF allows an 
implementation without any change to the model; the 
model only needs to be capable of exporting its state and 
restarting from the state modified by the EnKF. Conver- 
gence of the ensemble covariance to the state covariance 
is guaranteed in the large ensemble limit (3), (4), but a 
good approximation may require hundreds of ensemble 
members (2) because the covariance between physically 
distant variables is small, yet the sample covariance for a 
small ensemble has many large long-range terms. Local- 
ization techniques (5)-(7) improve the approximation by 
suppressing the long-range terms. 

In (8|-|T0|, we have proposed an alternative approach, 
the fast Fourier transform (FFT) EnKF. The FFT EnKF as- 
sumes that the state is a random field that is approximately 
homogeneous in space. Then the covariance matrix in the 
frequency domain can be well represented by its diagonal, 
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which gives a good approximation even for very small 
ensembles. However, the covariance is not represented 
well when it varies with location. The sample covariance 
can be used for the cross-covariance between different 
physical fields in the state [10], which, however, may 
again cause spurious long-range correlations. 

In this paper we extend the spectral approach to the 
wavelet EnKF, resulting in an automatic localization that 
varies in space adaptively. We also introduce a new tech- 
nique for automatic localization of the cross-covariances, 
based on a projection and a diagonalization in the spectral 
space. The efficiency of the new methods is demonstrated 
on an example. 

II. Related work 

Diagonal approximation of the covariance in the fre- 
quency space was proposed for weather fields (TTJ. 
Wavelets are well suited for approximation of meteoro- 
logical fields fT2| , and the diagonal approximation was 
extended to wavelet spaces |T3|-p5|. The Fourier domain 
KF (T6| is the KF applied to independent frequency 
modes. The Laplace operator represented by a diagonal 
matrix in the frequency space was used for a fast OSI 
(TO). The inverse of the Laplace operator was proposed 
as a covariance model [17], but higher negative powers 
(TO) yield better distributions. 

III. The ensemble Kalman filter 

The modeled quantity is the probability distribution of 
the state u, represented by an ensemble of simulation 
states {u\, . . . , un}. A data vector d is linked with the 
state u by the observation matrix H such that if the model 
and the data are correct, then d = Hu. The data error is 
assumed to have normal distribution with zero mean and 
known covariance R. When the data arrives, the ensemble 
is updated by 

u% = u k + Q N H T (HQ n H t + R)~\d + e k - Hu k ) , 

(1) 

where Qn is the sample covariance computed from 
the ensemble, and the data perturbation vectors e k are 
sampled from the data error distribution. The ensemble 
members are then advanced in time by the model until the 
next data vector is to be assimilated. See [2] for details. 



the entry of the vector u, corresponding to node X{. If the 
random field u is stationary, then the covariance matrix 



satisfies Q (x{ 



c (xi — Xj ) for some covariance 



function c, and matrix-vector multiplication v = Qu is 
the convolution, 



Figure. 1. Wavelet transform matrix with n — 64 and 5 
octaves, for the Coiflet 2 wavelet. 



IV. Orthogonal wavelet transform 
For vector u, denote u — Fu the transform 

u = [ue\t =1 , ug = fi-u, 

where fg are the rows of F and F is orthonormal, F~ x = 
F T . Matrices are then transformed by 

M = FMF~ 1 = FMF T . 

In the FFT EnKF, we use the Fourier sine transform. 
In the wavelet EnKF, we use the orthogonal wavelet 
transform, where the rows of F are fg = f( k ,j)> given by 
values of the family of functions (j)(2 k x — j) at the points 
x = i/n for a total n of composite indices £ = (k,j) 
such that < j < 2 k — 1. Each range of the indices 
(k,j) with a fixed k is called an octave. The function 
(j) is called the mother wavelet and is chosen so that the 
rows fi are orthonormal. The indices need to span full 
octaves, which restricts the dimension n to a power of 2 
(cf. Fig. [I]), though generalizations are possible. See (12), 
p8| for further details. We use WaveLab [19 ] to perform 
the wavelet transform u = Fu in Matlab. The complexity 
of the fast wavelet transform is only O (n), compared to 
O (n log n) for the FFT. 

In more than ID, the spectral transformations are ap- 
plied in each dimension separately, i.e., the basis functions 
are taken to be tensor products of ID basis functions. 
Unlike in the case of the FFT, the tensor product of 
wavelets creates a certain bias for coordinate directions; 
the impact of this bias on specific applications must be 
examined [13]. However, 2D wavelets do not seem to be 
practical yet. 

V. Spectral approximation of the covariance 
A. Single variable 

Consider first the case when the model state consists of 
one variable, in ID only. Denote by u(xi), i = 1, . . . , n 



v ( X i) = ^jQ X j) U ( X j) = ^2 U ( X j) C ( X i ~ X j) ' 

In the spectral domain, convolution becomes entry-by- 
entry multiplication of vectors, that is, the multiplication 
by a diagonal matrix, at least approximately. 

Let uik be the entries of the column vector Uk = Fuk, 
i.e., ensemble member k in the spectral space. Then we 
have the sample covariance in the spectral domain, 
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where u = -^ J2 k=1 Uk is the sample mean. Assuming 
that the covariance is primarily a function of the physical 
distance between the nodes Xi, we approximate the fore- 
cast covariance in the spectral domain by the diagonal 

D(u, u) of C(u, u), 
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(3) 



for i ^ j. 



B. Multiple variables 

When the model state consists of more than one 
variable, the covariance is split into blocks of cross- 
co variances between each variable. The diagonal blocks 
of the covariance can be approximated as in Off- 
diagonal blocks cannot, in general, be treated the same 
way because the meshes over which different variables 
are defined may not coincide. In this case, we define an 
interpolation operator P uv that projects a variable u to 
the mesh of variable v. While this matrix is rectangu- 
lar in general, we require an approximate left inverse, 
such as the Moore-Penrose pseudoinverse, P^ v , where 
P^ V P UV w I. In practice, more sophisticated and efficient 
methods with a sparse representation of the approximate 
left inverse are possible. In Section VII we consider 



only cross-covariances between identical grids, leaving a 
more detailed examination of these techniques for future 
research. 

By projecting u to the grid of v, it is possible to proceed 
as in the single variable case and construct an approximate 
cross-covariance, 

C(u, v) « PlC(P uv u, v) « PlF^d (P uv u, v) F. 



VI. Spectral EnKF 

A. Single variable 

First assume that the observation function H = I, that 
is, the whole state u is observed. The evaluation of the 
EnKF formula ([T]) in the frequency domain, with the di- 
agonal spectral approximation D(u,u) of the co variance 
becomes 

u a = u + D(u,u) (d(u,u) + p) 1 (3+e-u). (4) 

The analysis ensemble is obtained by the inverse trans- 
form at the end, u% = P _1 S^. Since D is diagonal, 
^ can be evaluated very efficiently in several important 
cases: (i) R is diagonal; then D + R is also diagonal and 
the evaluation of ^ reduces to term-by-term operations 
on vectors; (ii) R is the perturbed data sample covariance 
and we can use the Sherman-Morrison- Woodbury for- 
mula, which only requires solving a system of dimension 
equal to the ensemble size N; (Hi) R is approximated by 
the diagonal part of the sample covariance in the spectral 
domain, just like the state covariance. 

B. Multiple variables 

Consider the state with multiple variables and the 
covariance and the observation matrix in the block form. 
Using the spectral covariance estimation in Sec. |VJ 

QH T = [C(^,P»] « P+P" 1 \D(PiUi,H jU j\ P, 



Variable 1 



Variable 2 



with Pi the interpolation operator from mesh i to the 
observation grid, P^ the block diagonal matrix made 
up of P\ , • • • , Pt , and F consisting likewise of the 
spectral transform in each variable. Similarly, we have 



P, where the blocks are 
R in {!]) can be inverted 



HQH T « P- 1 D(HiU, HjU 

diagonal, so the term HQH T 
easily in the spectral domain. Since the cross covariances 
in the term QH T in {]]) are not involved in a matrix 
inversion, we can use another approximation, such the 
usual sample covariance 1 10 ], which may be more suitable 
when the field depends on the observed field non-locally, 
such as by advection. 

VII. Computational example 

We conclude with a simple example highlighting the 
advantages of the proposed method. We create a synthetic 
two-variable model state in one spatial dimension. Both 
variables are discretized on the same mesh of size 128 
over the domain [0,1]. The first variable is a simple 
Gaussian shape with random center, width, and height 
defined by 



u\(x) = /iexp(— (x — c) /w 2 ), 



(5) 



where c ~ A/"(0.3, 0.1 2 ), w ~ A/"(0.1, 0.01 2 ), and 
h ~ A/*(l,0.1 2 ). The second variable is made up of the 
sum of two components. The first is a smooth random 
field made up of a sum of sine functions with random 
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(a) (b) 
Figure. 2. A single realization of the stochastic model. 



amplitudes. The amplitudes are sampled from a Gaussian 
distribution with variance decaying as the inverse square 
of the frequency. The second component is identical to 
the first variable scaled by 0.3. The variables are chosen 
in this way so that the diagonal of Cov(ui,ui) is non- 
uniform with a peak at x = 0.3. Further, it is not obvious 
looking at a single realization (Figure [2]) that the two 
variables are correlated; however, it is easily verified that 
Cov(ui,U2) = 0.3 Cov(i£i, U\). 

Using a relatively small ensemble of size 10, one can 
see from Figure [3] that both the FFT and wavelet covari- 
ance estimates decay to zero with distance as expected. 
In comparison to a large sample covariance (size 1000), 
the spectral estimates offer a far better estimate of the 
true covariance than the traditional sample covariance 
when computed with the same ensemble. In addition, the 
wavelet covariance estimate correctly captures the spatial 
variability of the first variable's distribution. This is in 
contrast to the FFT estimate that smooths the peak located 
at 0.3 uniformly across the domain. 

In order to test the effectiveness of the spectral EnKF 
itself, we simulate an observation of the first variable with 
an observation grid that corresponds to the discretization 
of the model variables. We choose an observation gen- 
erated from §5§ with c = 0.4, w = 0.12, and h = 1.5. 



We apply the spectral EnKF described in Section [VI] with 
spectral transformations constructed from the discrete sine 
and Coiflet 2 wavelet transforms. We use same ensemble 
of size 10 from Figure [3] for each assimilation test and 
choose a small data covariance R = 0.01 2 / in order to 
force the assimilation to react visibly to the observation. 

Figure [4] shows the forecast and analysis of the first 
ensemble member for each method. The data is displayed 
in each figure to help gauge the accuracy of the response 
of the observed peak at x = 0.4. The traditional EnKF 
produces spurious noise near the peak of the forecast 
variable one, while the innovation in variable two is large 
throughout the domain rather than local to the observed 
peak. The FFT EnKF seems to react well to the observa- 
tion in the second variable, but the first variable exhibits 
spurious noise similar to a Gibb's effect throughout the 
domain. Finally, wavelet EnKF appears to provide the best 
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Figure. 3. In (b)-(d), covariance estimates with an ensemble 
size 10 are compared with the sample covariance of a large 
ensemble of size 1000 in (a). The color scale is the same in 
each figure with positive and negative correlations displayed in 
red and blue, respectively. 
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Figure. 4. The first (a-c) and second (d-f) variables of a single 
ensemble member for each algorithm compared. The analysis 
is displayed as a solid blue line along with the forecast as a 
dashed black line and observation as a dotted red line. 



analysis with very little spurious noise and a seemingly 
appropriate reaction to the observation. 

VIII. Conclusion 

Preliminary results indicate that the wavelet EnKF 
offers a vast improvement over the traditional assimilation 
techniques with similar ensemble sizes and without any 
localization needed. It requires only spatial information 
of the computational mesh, with no expert knowledge 
necessary to construct background covariances. Future 
work will analyze the effects of varying computational 



meshes as well as more complex observation functions 
on a range of operational computational models. 
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